% *** AutoPhase ***
% Created Date : 26 Oct 2013
% Modified Date : 26 Oct 2013
% Author : Sungjin Kim
% Description : Automatically correction of the phase from multi-coil
%               acquisition
% Arguments : 2D array with channels (data points by channel), ... 
%             c_args{1} - real(1), imaginary(2), or magnitude(3) part selection
%             c_args{2} - bandwidth of fid  
% ************************************************************************

%% BEGIN OF FUNCTION %%
function arr_phased_spectrum = autoPhase(arr_cmplxfid_data, c_args)
    
    interval_radian = 0.1;
    interval_freq = 1;
    size_data = size(arr_cmplxfid_data);
    nData = size_data(1);
    nChannel = size_data(2);
        
    % get SNR information
    arr_SNR = zeros([nChannel, 1]);
    for indx = 1:1:nChannel
        arr_spectra_data = fftshift(fft(arr_cmplxfid_data(:, indx)));
        arr_SNR(indx) = mrsSNR(arr_spectra_data, [1 501], 1, c_args);
    end
    
    [ref_SNR, ref_channel] = max(arr_SNR);
    for indx = 1:1:nChannel
        fprintf('The SNR of Channel %d : %f \n', indx, arr_SNR(indx));
    end
    fprintf('Reference coil number %d, Max SNR %f \n', ref_channel, ref_SNR);
    
    % get correlration value between reference data and others
    ref_channelData = arr_cmplxfid_data(:, ref_channel);
    ref_channelData = fftshift(fft(ref_channelData));
    arr_corrByzero = cell([nChannel,1]);  
    zero_correct = zeros([nChannel, 1]);

    % temporary value
    first_offsetFreq = 0;
    first_order_correction = 0;
    
    % 1 pass
    for indx_channel = 1:1:nChannel
            count_indx = 0;
            
            for indx_zero_order = 0:interval_radian:2*pi
                count_indx = count_indx + 1;
                arr_cmplxtemp = phaseShift(arr_cmplxfid_data(:, indx_channel), c_args{2}, ...
                                                indx_zero_order, first_order_correction, ...
                                                first_offsetFreq);
                temp = corrcoef(ref_channelData, arr_cmplxtemp); 
                arr_corrByzero{indx_channel}(count_indx) = real(temp(2,1));                                           
            end
            
            [value, indxVal] = max(arr_corrByzero{indx_channel});
            zero_correct(indx_channel) = (indxVal-1)*interval_radian;
            fprintf('The zero correlation value of channel %d is %f at %f radian\n', indx_channel, value, (indxVal-1)*interval_radian);

            arr_phased_spectrum(:, indx_channel) = phaseShift(arr_cmplxfid_data(:, indx_channel), c_args{2}, ...
                                    indxVal, first_order_correction, ...
                                    first_offsetFreq);        
    end
    
    % 2 pass
    for indx_channel = 1:1:nChannel
            count_indx = 0;
            
            for first_order_correction = 0:interval_radian:2*pi
                count_indx = count_indx + 1;
                arr_cmplxtemp = phaseShift(arr_cmplxfid_data(:, indx_channel), c_args{2}, ...
                                                zero_correct(indx_channel), first_order_correction, ...
                                                first_offsetFreq);
                temp = corrcoef(ref_channelData, arr_cmplxtemp); 
                arr_corrByzero{indx_channel}(count_indx) = real(temp(2,1));                                           
            end
            
            [value, indxVal] = max(arr_corrByzero{indx_channel});
            first_correct(indx_channel) = (indxVal-1)*interval_radian;
            fprintf('The first correlation value of channel %d is %f at %f radian\n', indx_channel, value, (indxVal-1)*interval_radian);

            arr_phased_spectrum(:, indx_channel) = phaseShift(arr_cmplxfid_data(:, indx_channel), c_args{2}, ...
                                    zero_correct(indx_channel), first_correct(indx_channel), ...
                                    first_offsetFreq);        
    end
    
%     % 3 pass
%     for indx_channel = 1:1:nChannel
%         count_indx = 0;
% 
%         for first_offsetFreq = 0:interval_freq:c_args{2}/4
%             first_offsetFreq = first_offsetFreq - c_args{2}/8;
%             count_indx = count_indx + 1;
%             arr_cmplxtemp = phaseShift(arr_cmplxfid_data(:, indx_channel), c_args{2}, ...
%                                             zero_correct(indx_channel), first_correct(indx_channel), ...
%                                             first_offsetFreq);
%             temp = corrcoef(ref_channelData, arr_cmplxtemp); 
%             arr_corrByzero{indx_channel}(count_indx) = real(temp(2,1));                                           
%         end
% 
%         [value, indxVal] = max(arr_corrByzero{indx_channel});
%         first_offset_correct(indx_channel) = (indxVal-1)*interval_freq;
%         fprintf('The offset correlation value of channel %d is %f at %f radian\n', indx_channel, value, (indxVal-1)*interval_freq);
% 
%         arr_phased_spectrum(:, indx_channel) = phaseShift(arr_cmplxfid_data(:, indx_channel), c_args{2}, ...
%                                 zero_correct(indx_channel), first_correct(indx_channel), ...
%                                 first_offset_correct(indx_channel));        
%     end
    
       
    return;
end % end of function